function meshes = BuildResampleMesh()
% This function computes a set of local meshes that can be used to
% resample a finite element space on an element. Typically, this is
% useful for the visualization.
%
% The returned variable meshes is a three index array as follows:
% first index:  number of spatial dimensions
% second index: global polynomial degree
% third index:  local polynomial degree
%
% The distinction between global and local polynomial degrees allows
% to represent (interpolate), for instance, a 2D quadratic triangle as
% a collection of three linear triangles. Notice that not all the
% combinations global/local polynomial degree make sense; however, at
% least the two cases
%   local = global
%   local = 1
% always make sense.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 1D                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 1              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 ];
t = [ 1 ;
      2 ];
PV_cell_type = 3;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,1,1) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 2              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0.0 , 0.5 , 1 ];
t = [ 1 , 2 ;
      2 , 3 ];
PV_cell_type = 3;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,2,1) = s;

t = [ 1 ;
      2 ;
      3 ];
PV_cell_type = 21;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,2,2) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 3              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0.0 , 1/3 , 2/3 , 1 ];
t = [ 1 , 2 , 3 ;
      2 , 3 , 4 ];
PV_cell_type = 3;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,3,1) = s;

t = [ 1 ;
      2 ;
      3 ;
      4 ];
PV_cell_type = -1;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,3,3) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 4              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0.0 , 0.25 , 0.5 , 0.75 , 1 ];
t = [ 1 , 2 , 3 , 4 ;
      2 , 3 , 4 , 5 ];
PV_cell_type = 3;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,4,1) = s;

t = [ 1 , 3 ;
      2 , 4 ;
      3 , 5 ];
PV_cell_type = 21;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,4,2) = s;

t = [ 1 ;
      2 ;
      3 ;
      4 ;
      5 ];
PV_cell_type = -1;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(1,4,4) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 2D                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 1              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 ;
      0 , 0 , 1 ];
t = [ 1 ;
      2 ;
      3 ];
PV_cell_type = 5;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,1,1) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 2              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 , 0.5 , 0.5 , 0.0 ;
      0 , 0 , 1 , 0.0 , 0.5 , 0.5 ];
t = [ 1 , 4 , 6 , 4 ;
      4 , 2 , 5 , 5 ;
      6 , 5 , 3 , 6 ];
PV_cell_type = 5;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,2,1) = s;

t = [ 1 ;
      2 ;
      3 ;
      4 ;
      5 ;
      6 ];
PV_cell_type = 22;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,2,2) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 3              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 , 1/3 , 2/3 , 2/3 , 1/3 , 0.0 , 0.0 , 1/3 ;
      0 , 0 , 1 , 0.0 , 0.0 , 1/3 , 2/3 , 2/3 , 1/3 , 1/3 ];
t = [ 1 , 5 , 8 ,  4 ,  4 ,  5 , 10 , 10 ,  9 ;
      4 , 2 , 7 , 10 ,  5 ,  6 ,  6 ,  7 , 10 ;
      9 , 6 , 3 ,  9 , 10 , 10 ,  7 ,  8 ,  8 ];
PV_cell_type = 5;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,3,1) = s;

t = [ 1 ;
      2 ;
      3 ;
      4 ;
      5 ;
      6 ;
      7 ;
      8 ;
      9 ;
     10 ];
PV_cell_type = -1;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,3,3) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 4              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 , 1/4 , 2/4 , 3/4 , 3/4 , 2/4 , 1/4 , 0.0 , 0.0 , 0.0 , 1/4 , 2/4 , 1/4 ;
      0 , 0 , 1 , 0.0 , 0.0 , 0.0 , 1/4 , 2/4 , 3/4 , 3/4 , 2/4 , 1/4 , 1/4 , 1/4 , 2/4 ];
t = [  1 ,  6 , 10 ,  4 ,  4 ,  5 ,  5 ,  6 , 14 , 14 , 15 , 15 , 11 , 13 , 12 , 13 ;
       4 ,  2 ,  9 , 13 ,  5 , 14 ,  6 ,  7 ,  7 ,  8 ,  8 ,  9 , 15 , 15 , 13 , 14 ;
      12 ,  7 ,  3 , 12 , 13 , 13 , 14 , 14 ,  8 , 15 ,  9 , 10 , 10 , 11 , 11 , 15 ];
PV_cell_type = 5;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,4,1) = s;

t = [  1 ,  5 , 11 ,  5 ;
       5 ,  2 ,  8 ,  8 ;
      11 ,  8 ,  3 , 11 ;
       4 ,  6 , 15 , 14 ;
      13 ,  7 ,  9 , 15 ;
      12 , 14 , 10 , 13 ];

PV_cell_type = 22;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,4,2) = s;

t = [ 1 ;
      2 ;
      3 ;
      4 ;
      5 ;
      6 ;
      7 ;
      8 ;
      9 ;
     10 ;
     11 ;
     12 ;
     13 ;
     14 ;
     15 ];
PV_cell_type = -1;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,4,4) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 5              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 , 1/5 , 2/5 , 3/5 , 4/5 , 4/5 , 3/5 , 2/5 , 1/5 , 0.0 , 0.0 , 0.0 , 0.0 , 1/5 , 2/5 , 3/5 , 2/5 , 1/5 , 1/5 ;
      0 , 0 , 1 , 0.0 , 0.0 , 0.0 , 0.0 , 1/5 , 2/5 , 3/5 , 4/5 , 4/5 , 3/5 , 2/5 , 1/5 , 1/5 , 1/5 , 1/5 , 2/5 , 3/5 , 2/5 ];
t = [  1 ,  4 ,  4 ,  4 ,  5 ,  6 ,  6 ,  6 ,  7 , 15 , 15 , 16 , 17 , 17 , 17 , 18  , 14 , 21 , 21 , 21 , 19 , 13 , 13 , 20 , 12 ;
       4 , 16 , 17 ,  5 ,  6 , 18 ,  8 ,  7 ,  2 , 16 , 21 , 17 , 19 , 18 ,  9 ,  8  , 21 , 20 , 19 , 10 ,  9 , 20 , 11 , 10 , 11 ;
      15 , 15 , 16 , 17 , 17 , 17 , 18 ,  8 ,  8 , 21 , 14 , 21 , 21 ,  9 , 19 ,  9  , 13 , 13 , 10 , 20 , 10 , 11 , 12 , 11 ,  3 ];
PV_cell_type = 5;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(2,5,1) = s;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 3D                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 1              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 , 0 ;
      0 , 0 , 1 , 0 ;
      0 , 0 , 0 , 1 ];
t = [ 1 ;
      2 ;
      3 ;
      4 ];
PV_cell_type = 10;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(3,1,1) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 2              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = [ 0 , 1 , 0 , 0 , 1/2 , 1/2 , 0.0 , 0.0 , 1/2 , 0.0 ;
      0 , 0 , 1 , 0 , 0.0 , 1/2 , 1/2 , 0.0 , 0.0 , 1/2 ;
      0 , 0 , 0 , 1 , 0.0 , 0.0 , 0.0 , 1/2 , 1/2 , 1/2 ];

% The main building block for the connectivity is now how to cut a
% regular prism into three tetrahedra. Assuming that the lower basis
% of the prism has vertices [1 2 3] and the upper one [4 5 6] (with 4
% laying "above" 1, 5 "above" 2 and 6 "above" 3), the prism can be cut
% in one of these two ways:
%
%      [ 1  2  3 ]          [ 1  2  3 ]
%  up: [ 2  4  6 ]    down: [ 2  4  6 ]
%      [ 3  5  4 ]          [ 3  5  4 ]
%      [ 4  6  2 ]          [ 4  3  5 ]

% An often needed object is a cube to which one tetrahedron has been
% removed: we call it CutCube. Assuming the numbering [ 1 2 3 4 ] for
% the lower face and [5 6 7] for the upper one (the vertex "above" 3
% is missing), we have
CutCube = [ 1 , 2 , 4 , 2 , 2 ;
            2 , 5 , 7 , 3 , 6 ;
            4 , 6 , 5 , 4 , 3 ;
            5 , 7 , 2 , 7 , 7 ];
% A complete cube can be obtained adding the missing tetrahedra (the
% two faces of the cube are now [1 2 3 4] and [5 6 7 8], so that the
% numbering of the upper face is different from the one used in
% CutCube and a renumbering is required).
verts = [1 2 3 4 5 6 8];
Cube = [ verts(CutCube) , [7;6;8;3] ];

% Fist cut three tetrahedra at the three external corners
t1 = [ 5 ,  7 ,  8 ;
       2 ,  6 ,  9 ;
       6 ,  3 , 10 ;
       9 , 10 ,  4 ];

% Then the remaining part is a CutCube with vertices [1 5 6 7 8 9 10]
verts = [1 5 6 7 8 9 10];
t2 = verts(CutCube);

t = [ t1 , t2 ];
PV_cell_type = 10;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(3,2,1) = s;

t = [ 1 ;
      2 ;
      3 ;
      4 ;
      5 ;
      6 ;
      7 ;
      8 ;
      9 ;
     10 ];
PV_cell_type = 24;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(3,2,2) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 3              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p = [ 0 , 1 , 0 , 0 , 1/3 , 2/3 , 2/3 , 1/3 , 0.0 , 0.0 , 1/3 , 0.0 , 2/3 , 0.0 , 1/3 , 1/3 , 0.0 , 0.0 , 1/3 , 0.0 ;
      0 , 0 , 1 , 0 , 0.0 , 0.0 , 1/3 , 2/3 , 2/3 , 1/3 , 1/3 , 0.0 , 0.0 , 2/3 , 0.0 , 1/3 , 1/3 , 0.0 , 0.0 , 1/3 ;
      0 , 0 , 0 , 1 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1/3 , 1/3 , 1/3 , 1/3 , 1/3 , 1/3 , 2/3 , 2/3 , 2/3 ];

% Fist cut three tetrahedra at the three external corners
t1 = [  6 ,  9 , 18 ;
        2 ,  8 , 19 ;
        7 ,  3 , 20 ;
       13 , 14 ,  4 ];

% Now cut three CutCube-s
verts1 = [ 5  6  7 11 15 13 16];
verts2 = [10 11  8  9 17 16 14];
verts3 = [12 15 16 17 18 19 20];
t2 = [ verts1(CutCube) , verts2(CutCube) , verts3(CutCube) ];

% Now three tetrahedra facing the x+y+z=1 face
t3 = [ 11 , 15 , 17 ;
        7 , 13 , 16 ;
        8 , 16 , 14 ;
       16 , 19 , 20 ];

% Finally a cube at the origin
verts = [1 5 11 10 12 15 16 17];
t4 = verts(Cube);

t = [ t1 , t2 , t3 , t4 ];
PV_cell_type = 10;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(3,3,1) = s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           degree 4              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p1 = [ 0 , 1 , 0 , 0 , 1/4 , 2/4 , 3/4 , 3/4 , 2/4 , 1/4 , 0.0 , 0.0 , 0.0 , 1/4 , 2/4 , 1/4 ;
       0 , 0 , 1 , 0 , 0.0 , 0.0 , 0.0 , 1/4 , 2/4 , 3/4 , 3/4 , 2/4 , 1/4 , 1/4 , 1/4 , 2/4 ;
       0 , 0 , 0 , 1 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ];
p2 = [ 0.0 , 3/4 , 0.0 , 1/4 , 2/4 , 2/4 , 1/4 , 0.0 , 0.0 , 1/4 ;
       0.0 , 0.0 , 3/4 , 0.0 , 0.0 , 1/4 , 2/4 , 2/4 , 1/4 , 1/4 ;
       1/4 , 1/4 , 1/4 , 1/4 , 1/4 , 1/4 , 1/4 , 1/4 , 1/4 , 1/4 ];
p3 = [ 0.0 , 2/4 , 0.0 , 1/4 , 1/4 , 0.0 ;
       0.0 , 0.0 , 2/4 , 0.0 , 1/4 , 1/4 ;
       2/4 , 2/4 , 2/4 , 2/4 , 2/4 , 2/4 ];
p4 = [ 0.0 , 1/4 , 0.0 ;
       0.0 , 0.0 , 1/4 ;
       3/4 , 3/4 , 3/4 ];
p = [ p1 , p2 , p3 , p4 ];

% Fist cut three tetrahedra at the three external corners
t1 = [  7 , 11 , 33 ;
        2 , 10 , 34 ;
        8 ,  3 , 35 ;
       18 , 19 ,  4 ];

%% Now cut six CutCube-s
verts1 = [ 6  7  8 15 21 18 22];
verts2 = [14 15  9 16 26 22 23];
verts3 = [12 16 10 11 24 23 19];
verts4 = [20 21 22 26 30 28 31];
verts5 = [25 26 23 24 32 31 29];
verts6 = [27 30 31 32 33 34 35];
t2 = [ verts1(CutCube) , verts2(CutCube) , verts3(CutCube) , verts4(CutCube) , verts5(CutCube) , verts6(CutCube) ];

% Now seven tetrahedra facing the x+y+z=1 face
t3 = [ 15 , 16 , 21 , 26 , 24 , 30 , 32 ;
        8 ,  9 , 18 , 22 , 23 , 28 , 31 ;
        9 , 10 , 22 , 23 , 19 , 31 , 29 ;
       22 , 23 , 28 , 31 , 29 , 34 , 35 ];

% Finally four cubes at the origin
verts1 = [ 1  5 14 13 17 20 26 25];
verts2 = [ 5  6 15 14 20 21 22 26];
verts3 = [13 14 16 12 25 26 23 24];
verts4 = [17 20 26 25 27 30 31 32];
t4 = [ verts1(Cube) , verts2(Cube) , verts3(Cube) , verts4(Cube) ];

t = [ t1 , t2 , t3 , t4 ];
PV_cell_type = 10;

s = struct('p',p,'t',t,'PV_cell_type',PV_cell_type);
meshes(3,4,1) = s;

return

